#install.packages("BiocManager")
#BiocManager::install("Seurat")
#BiocManager::install("ggplot2")
#BiocManager::install("sctransform")
library(Seurat)
library(ggplot2)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ape)
library(cowplot)
library(Matrix)
library(EnhancedVolcano)
## Loading required package: ggrepel
#find out where you are
getwd()
## [1] "/Users/whippoorwill/Desktop/Code/scrnaseq/Rmd"
#Specify where your matrix files are
dir= "/Users/whippoorwill/Desktop/Sequencing/LD_AVM02/"
datafolder = "Data/Seurat"
filename = "Microglia_BC_Macrophages_subset.RData"
organism = "Mouse"
defile = "Macrophage_only_all_markers.csv"
m = c("nCount_RNA","nFeature_RNA","percent.mito","percent.ribo")
Plotfolder = "Plots"
if(organism == "Mouse"){library(org.Mm.eg.db)}
## Loading required package: AnnotationDbi
## Loading required package: stats4
## Loading required package: BiocGenerics
## Loading required package: parallel
##
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB
## The following object is masked from 'package:Matrix':
##
## which
## The following objects are masked from 'package:dplyr':
##
## combine, intersect, setdiff, union
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, basename, cbind, colnames,
## dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
## grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
## rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
## union, unique, unsplit, which, which.max, which.min
## Loading required package: Biobase
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
## Loading required package: IRanges
## Loading required package: S4Vectors
##
## Attaching package: 'S4Vectors'
## The following object is masked from 'package:Matrix':
##
## expand
## The following objects are masked from 'package:dplyr':
##
## first, rename
## The following object is masked from 'package:base':
##
## expand.grid
##
## Attaching package: 'IRanges'
## The following objects are masked from 'package:dplyr':
##
## collapse, desc, slice
##
## Attaching package: 'AnnotationDbi'
## The following object is masked from 'package:dplyr':
##
## select
##
if(organism == "Human"){library(org.Hs.eg.db)}
if(organism == "Zebrafish"){library(org.Dr.eg.db)}
Load in your filtered dataset
load(file.path(dir,datafolder,filename))
Other differential expression methods:
base = '0' #default cluster
column = "seurat_clusters" #can pick any metadata column
Idents(sobject) = column
clusters = levels(sobject@meta.data[,column])
clusters = clusters[clusters != base]
for (cluster in clusters){
markers_all <- FindMarkers(
object = sobject,
ident.1 = cluster,
ident.2 = base,
only.pos = FALSE,
min.pct = 0.10, #gene must be present in 10% of the cells in the cluster
logfc.threshold = 0,
test.use = "MAST")
dim(markers_all)
head(markers_all)
write.csv(markers_all,file = file.path(dir,"Spreadsheets",paste0("markers",cluster,"_vs_",base,".csv")))
}
Or you can run “bulk-seq” like analyses based on your original sample IDs:
#can pick any metadata column
column = "sample_description"
#default cluster
cluster1 = 'Deprived-P5'
#cluster of interest
cluster2 = "Control-P5"
Idents(sobject) = column
markers = FindMarkers(sobject,
ident.1=cluster1,
ident.2 = cluster2,
only.pos=F,
logfc.threshold = 0.0,
min.pct = 0.1,
test.use = "MAST")
write.csv(markers,file.path(dir,"Spreadsheets",paste0("markers_",cluster1,"_vs_",cluster2,".csv")))
You can also make plots with only subsets of cells: Best for making umap plots where both conditions/individuals have equal representation.
#first check how many cells are in each group so you don't pick a number more than the min
column = "sample_description"
ncells = 2000
genes = c("Cx3cr1","P2ry12","Spp1","Nrxn2","Ifitm3","Mki67","Pf4")
features = m
categories = c("sample_description","age","condition","celltypecluster")
name = "macrophage_equalcells"
Function to print multiple graphs:
PrintSeuratGraph = function(namecard = "a",seurat_object = sobject,graphtype = "feature",feature = NULL,group = NULL,split=NULL,cellnames=NULL){
if (!is.null(cellnames)){
Idents(seurat_object) = cellnames[1]
cells = colnames(seurat_object)[Idents(seurat_object) %in% cellnames[2:length(cellnames)]]}
else {cells = cellnames}
if (graphtype == "feature"){
graph = FeaturePlot(seurat_object,features = feature,split.by = split, cells = cells,cols = c("lightyellow","darkred"))
}
if (graphtype == "violin"){
graph = VlnPlot(seurat_object,features = feature, pt.size = 0.1, idents = cellnames[2:length(cellnames)],group.by = group, split.by = split)
}
if (graphtype == "dim"){
graph = DimPlot(seurat_object,cells = cells, group.by = group, split.by = split)
}
name = paste0(feature,"_",graphtype,namecard,".eps")
graph
setEPS()
postscript(file.path(dir,Plotfolder,name))
print(graph)
dev.off()
}
table(sobject$sample_description)
##
## Control-P5 Control-P7 Deprived-P5 Deprived-P7
## 2975 4125 2566 2743
cellnames = sobject@meta.data[,column]
names(cellnames) = colnames(sobject)
groups = levels(cellnames)
newcellnames = NULL
for (group in groups){
cells = sample((cellnames)[cellnames == group],ncells)
newcellnames = c(newcellnames,cells)
}
#feature plots
for(feature in c(genes,features)){
PrintSeuratGraph(namecard = name,graphtype = "feature",feature = feature,cellnames = newcellnames)
}
#split feature plots by individual
for(feature in c(features)){
PrintSeuratGraph(namecard = paste0(name,"_split"),graphtype = "feature",feature = feature,split = "sample_description",cellnames = newcellnames)
}
#dim plots for clustering
for(group in categories){
PrintSeuratGraph(namecard = name,graphtype = "dim",group = group, feature = group,cellnames = newcellnames)
}
#violin plots
for(feature in c(genes,features)){
PrintSeuratGraph(namecard = name,graphtype = "violin",feature = feature,group = "seurat_clusters",cellnames = newcellnames)
}
You can save umap and pc embeddings to load into other programs like Loupe and anything in Python like scvelo (velocity)
#Edit name only if one project has multiple embeddings
name = paste0(Project(sobject),"_",name)
#Save umap embeddings
umap =as.data.frame(sobject@reductions$umap@cell.embeddings)
umap$Cellname = rownames(umap)
umap = umap[,order(colnames(umap))]
write.csv(umap,file = file.path(dir,"Annotation",paste0(name,"_umap_embed.csv")),row.names = F)
#save pcs
pc =as.data.frame(sobject@reductions$pc@cell.embeddings)
pc$Cellname = rownames(pc)
pc = pc[,order(colnames(pc))]
write.csv(pc,file = file.path(dir,"Annotation",paste0(name,"_pc_embed.csv")),row.names = F)
#save metadata
meta = sobject@meta.data
meta$Cellname = rownames(meta)
meta = meta[,c(ncol(meta),1:(ncol(meta)-1))]
write.csv(pc,file = file.path(dir,"Annotation",paste0(name,"_metadata.csv")),row.names = F)
#For PanoView (in Python, a way to mathematically validate # of clusters found)
annotation = sobject$sample_description #or any column
x = as.matrix(GetAssayData(sobject,slot = "counts"))
x = x[rownames(x) %in% VariableFeatures(sobject),]
write.csv(x,file = file.path(dir,"Annotation",paste0(name,"_counts.csv")))
write.csv(annotation,file = file.path(dir,"Annotation",paste0(name,"_annotation.csv")))
You can add annotations from the annotationDB database:
de = read.csv(file.path(dir,"Spreadsheets",defile),stringsAsFactors = F) #any spreadsheet with gene symbols or other identifiers
if (organism == "Mouse"){db = org.Mm.eg.db}
if (organism == "Human"){db = org.Hs.eg.db}
if (organism == "Zebrafish"){db = org.Dr.eg.db}
ids=de$gene
fromKey="SYMBOL" #must match the ids - could also be ensembl ID
toKey=c("GENENAME","ENSEMBL","UNIPROT") #whatever annotation you want to add - find with keytypes(db)
selRes<-AnnotationDbi::select(db,keys=ids,keytype=fromKey,columns=c(fromKey,toKey))
## 'select()' returned many:many mapping between keys and columns
x=selRes[match(ids,selRes[,1]),1:(length(toKey)+1)]
identical(x$SYMBOL,de$gene)
## [1] TRUE
de$GeneName = x$GENENAME
de$Ensembl = x$ENSEMBL
de$Uniprot = x$UNIPROT
Volcano Plot
Set your parameters
#Minimum fold change (i.e. 1.15 = 15% increase)
minfc = 1.15
#Max adj. p value
alpha = 1e-25
#Clusters selected
categories = levels(as.factor(de$cluster))
#Genes to highlight
ngenes = 20
Set up the spreadsheet correctly
colnames(de)[8] = "Gene"
newlist = list()
clusters = levels(as.factor(de$cluster))
#Split by cluster
i = 1
for (cluster in clusters){
newlist[[cluster]] = de[de$cluster == cluster,]
i = i+1
}
#select a single cluster
for (category in categories){
fc = newlist[[category]]
fc = fc[!is.na(fc$avg_logFC),]
colorkeysdown = fc$Gene[fc$avg_logFC < -log2(minfc) & fc$p_val_adj < alpha]
colorkeysup = fc$Gene[fc$avg_logFC > log2(minfc) & fc$p_val_adj < alpha]
#Either highlight specific genes or pick the top genes in colorkeysup/down
top = fc[fc$p_val_adj<alpha,]
top = top[order(top$avg_logFC),"Gene"]
highlight = c(head(top,ngenes),tail(top,ngenes))
allcolors = rep("darkgrey",length(fc$Gene))
names(allcolors) = fc$Gene
allcolors[names(allcolors) %in% colorkeysdown] = "blue"
allcolors[names(allcolors) %in% colorkeysup]= "red"
allcolors[names(allcolors) %in% highlight]= "yellow"
names(allcolors)[allcolors == "yellow"] = "labelled"
names(allcolors)[allcolors == "red"] = "u"
names(allcolors)[allcolors == "darkgrey"] = "-"
names(allcolors)[allcolors == "blue"] = "d"
setEPS()
postscript(file.path(dir,"Plots",paste0("Volcano_",category,".eps")))
print(EnhancedVolcano(fc,
lab = fc$Gene,
x = 'avg_logFC',
y = 'p_val_adj',
xlim = c(-3, 3),
title = category,
subtitle = "",
drawConnectors = F,
legendPosition = 'right',
legendVisible = F,
pCutoff = alpha,
FCcutoff = log2(minfc),
selectLab = highlight,
transcriptPointSize = 1.5,
transcriptLabSize = 2.0,
col=c('black', 'black', 'black', 'red3'),
colCustom = allcolors,
gridlines.major = F,
gridlines.minor = F,
colAlpha = 1))
dev.off()
print(EnhancedVolcano(fc,
lab = fc$Gene,
x = 'avg_logFC',
y = 'p_val_adj',
xlim = c(-3, 3),
title = category,
subtitle = "",
drawConnectors = F,
legendPosition = 'right',
legendVisible = F,
pCutoff = alpha,
FCcutoff = log2(minfc),
selectLab = highlight,
transcriptPointSize = 1.5,
transcriptLabSize = 2.0,
col=c('black', 'black', 'black', 'red3'),
colCustom = allcolors,
gridlines.major = F,
gridlines.minor = F,
colAlpha = 1))
}
## Warning: One or more P values is 0. Converting to minimum possible value...
## Warning: One or more P values is 0. Converting to minimum possible value...
## Warning: One or more P values is 0. Converting to minimum possible value...
## Warning: One or more P values is 0. Converting to minimum possible value...
## Warning: One or more P values is 0. Converting to minimum possible value...
## Warning: One or more P values is 0. Converting to minimum possible value...
## Warning: One or more P values is 0. Converting to minimum possible value...
## Warning: One or more P values is 0. Converting to minimum possible value...
## Warning: One or more P values is 0. Converting to minimum possible value...
## Warning: One or more P values is 0. Converting to minimum possible value...
## Warning: One or more P values is 0. Converting to minimum possible value...
## Warning: One or more P values is 0. Converting to minimum possible value...
## Warning: One or more P values is 0. Converting to minimum possible value...
## Warning: One or more P values is 0. Converting to minimum possible value...
You can find what phase of the cell cycle each cell is in:
Cell cycle regression specifically for mouse genes right now - can adjust slightly for human/zebrafish data by changing the capitalization of the spreadsheet
cc.genes <- readLines(con = "../Spreadsheets/regev_lab_cell_cycle_genes.txt")
genesofinterest = c("Pcna","Top2a","Mcm6","Mki67")
#for zebrafish and mice, change gene names to lowercase
if (organism %in% c("Zebrafish","Mouse")){cc.genes<-tolower(cc.genes)}
#for mice, capitalize the first letter of each gene
if (organism == "Mouse"){
cc.genes<-unname(sapply(cc.genes,function(x){
x<-paste0(toupper(substr(x,start = 1,stop = 1)),substring(x,first=2))
}))
}
#assign each gene to "S" phase or "G2M" phase
s.genes <- cc.genes[1:43]
g2m.genes <- cc.genes[44:97]
#score each cell by gene expression on this subset of genes
sobject <- CellCycleScoring(object = sobject, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)
## Warning: The following features are not present in the object: Mlf1ip, not
## searching for symbol synonyms
## Warning: The following features are not present in the object: Fam64a, Hn1, not
## searching for symbol synonyms
# view cell cycle scores and phase assignments
head(x = sobject@meta.data)
## orig.ident nCount_RNA nFeature_RNA percent.mito
## AAACCCAAGTCTACCA-LD5LD LD5LD 22762 4959 2.455847
## AAACGAAAGCGTGTCC-LD5LD LD5LD 10098 3218 3.069915
## AAACGAAAGGTCTTTG-LD5LD LD5LD 13770 4173 2.222222
## AAACGAATCGGTAGAG-LD5LD LD5LD 19687 4451 5.338548
## AAACGCTAGGTACATA-LD5LD LD5LD 13015 3875 3.234729
## AAACGCTAGTAGCATA-LD5LD LD5LD 12825 3816 4.421053
## percent.ribo sample_description age condition sex
## AAACCCAAGTCTACCA-LD5LD 19.299710 Deprived-P5 P5 Deprived Female
## AAACGAAAGCGTGTCC-LD5LD 22.291543 Deprived-P5 P5 Deprived Male
## AAACGAAAGGTCTTTG-LD5LD 12.541757 Deprived-P5 P5 Deprived Female
## AAACGAATCGGTAGAG-LD5LD 20.150353 Deprived-P5 P5 Deprived Female
## AAACGCTAGGTACATA-LD5LD 9.880907 Deprived-P5 P5 Deprived Female
## AAACGCTAGTAGCATA-LD5LD 19.555556 Deprived-P5 P5 Deprived Female
## nCount_SCT nFeature_SCT SCT_snn_res.0.5 SCT_snn_res.1.5
## AAACCCAAGTCTACCA-LD5LD 15123 4774 4 10
## AAACGAAAGCGTGTCC-LD5LD 13442 3223 0 0
## AAACGAAAGGTCTTTG-LD5LD 13974 4172 3 2
## AAACGAATCGGTAGAG-LD5LD 15110 4409 4 10
## AAACGCTAGGTACATA-LD5LD 13714 3874 2 4
## AAACGCTAGTAGCATA-LD5LD 13672 3816 0 0
## SCT_snn_res.1 seurat_clusters celltype celltypecluster
## AAACCCAAGTCTACCA-LD5LD 10 0 Microglia Microglia-0
## AAACGAAAGCGTGTCC-LD5LD 0 0 Microglia Microglia-0
## AAACGAAAGGTCTTTG-LD5LD 1 0 Microglia Microglia-0
## AAACGAATCGGTAGAG-LD5LD 10 0 Microglia Microglia-0
## AAACGCTAGGTACATA-LD5LD 8 2 Microglia Microglia-2
## AAACGCTAGTAGCATA-LD5LD 0 0 Microglia Microglia-0
## SCT_snn_res.0.1 S.Score G2M.Score Phase old.ident
## AAACCCAAGTCTACCA-LD5LD 0 -0.3239106 -0.4315991 G1 G1
## AAACGAAAGCGTGTCC-LD5LD 0 -0.2204614 -0.3003426 G1 G1
## AAACGAAAGGTCTTTG-LD5LD 0 0.3200799 -0.3523386 S S
## AAACGAATCGGTAGAG-LD5LD 0 -0.2900263 -0.4606724 G1 G1
## AAACGCTAGGTACATA-LD5LD 2 -0.1918543 -0.3783492 G1 G1
## AAACGCTAGTAGCATA-LD5LD 0 -0.1110625 -0.3937208 G1 G1
RidgePlot(object = sobject, features = genesofinterest)
## Picking joint bandwidth of 0.0685
## Picking joint bandwidth of 0.123
## Picking joint bandwidth of 0.0844
## Picking joint bandwidth of 0.135
FeaturePlot(sobject,"S.Score")
FeaturePlot(sobject,"G2M.Score")
#save phase plot
setEPS()
postscript(file.path(dir,"Plots",paste0(name,"_phasechart.eps")))
DimPlot(sobject,group.by = "Phase")
dev.off()
## quartz_off_screen
## 2
#save seurat object with cell cycle scoring
save(sobject,file = file.path(dir,datafolder,filename)))
You can regress out the phase genes
regress = c("S.Score","G2M.Score","nCount_RNA","percent.mito")
columns = c("sample_description","seurat_clusters","celltypecluster","condition","age","Phase")
#run PCA on only the cell cycle genes
sobject2 <- RunPCA(object = sobject, pc.genes = c(s.genes, g2m.genes), do.print = FALSE, maxit=10000)
#cell cycle only
PCAPlot(object = sobject2)
#original
PCAPlot(object = sobject)
#regress out cell cycle variables
sobject2 <- ScaleData(object = sobject2, features = VariableFeatures(sobject), vars.to.regress = regress)
# Now, a PCA on the variable genes no longer returns components associated
# with cell cycle
sobject2 <- RunPCA(object = sobject2, features = VariableFeatures(sobject), genes.print = 10)
PCAPlot(object = sobject2)
#remake umap with new calculations without cell cycle
sobject2@reductions$UMAP<-NULL
sobject2<-RunUMAP(sobject2,reduction = "pca",dims = 1:30, verbose = F)
sobject2<-FindNeighbors(sobject2,dims=1:30,verbose=F)
set.seed(1)
sobject2<-FindClusters(sobject2,verbose=F,resolution = 0.5)
load(file.path(dir,datafolder, paste0(name,"_cellcycleregressed.RData")))
#make plots without cell cycle
for (column in columns){
print(DimPlot(object = sobject2, group.by=column, pt.size=0.5,label = T))
}
## Warning: Using `as.character()` on a quosure is deprecated as of rlang 0.3.0.
## Please use `as_label()` or `as_name()` instead.
## This warning is displayed once per session.
#save seurat object with cell cycle scoring
save(sobject2,file = file.path(dir,datafolder,paste0(name,"_cellcycleregressed.RData")))
Barplot for any two (or more) categories
#Pick metadata columns
clustercolumn = "seurat_clusters"
samplecolumn = "sample_description"
#pick a reasonable number of cells per sample to normalize by
ncells = 2000
cols = c("blue","red")
#If you want to only compare particular samples/conditions, split further by another metadata column:
split = "age"
splitby = levels(as.factor(sobject[[split]][,1]))
#Make a table and normalize
r = table(sobject[[clustercolumn]][,1],sobject[[samplecolumn]][,1])
#Split the table
for (i in splitby){
x = grep(i,colnames(r))
t = r[,x]
#remove any clusters that don't have cells
t = t[rowSums(t)>0,]
#normalize by sample
t = apply(t,MARGIN = 2,function(x)x/sum(x))
t = round(t*ncells,0)
#convert to percents for each cluster
t = apply(t,MARGIN = 1,function(x)x/sum(x))
t = round(t*100,2)
setEPS()
postscript(file.path(dir,"Plots",paste0(name,i,"barplot.eps")))
barplot(t, main="Cluster composition by percent of celltype",
xlab="Cluster", ylab = "% of cluster", ylim = c(0,100), col=cols,axisnames = T,
width = .2,xlim = c(0,5),legend = rownames(t), space = 0.6,cex.names = 0.8,axis.lty = 1)
dev.off()
print(barplot(t, main="Cluster composition by percent of celltype",
xlab="Cluster", ylab = "% of cluster", ylim = c(0,100), col=cols,axisnames = T,
width = .2,xlim = c(0,5),legend = rownames(t), space = 0.6,cex.names = 0.8,axis.lty = 1))
}
## [1] 0.22 0.54 0.86 1.18 1.50 1.82 2.14
## [1] 0.22 0.54 0.86 1.18 1.50 1.82 2.14
Contingency tables for barplots. Note that this makes more sense if done with only 2 samples, as only one statistic will be calculated per cluster.
clustercolumn = "seurat_clusters"
samplecolumn = c("sample_description")
all = table(sobject[[clustercolumn]][,1],sobject[[samplecolumn]][,1])
array = as.data.frame(cbind("Cluster" = 1, "Sample1" = 1, "Sample2" = 1, "p-value" = 1, "Cramer's V" =1))
for (i in 1:(ncol(all)-1)){
for (k in (i+1):ncol(all)){
partial = all[,c(i,k)]
total = colSums(partial)
clusters = rownames(partial)
test = list()
chi = chisq.test(partial)
for (cluster in clusters){
x = rbind(partial[cluster,],total-partial[cluster,])
rownames(x) = c(cluster,"total")
chi = chisq.test(x)
f = sqrt(chi$statistic/sum(x))
test[[cluster]] = chi
test[[cluster]]$f = f
pval = test[[cluster]]$p.value
Cramer = sqrt(test[[cluster]]$statistic/sum(all[cluster,]))
out = c(cluster,colnames(partial),pval,Cramer)
#print out the p-value and cramer's v (closer to 1 = higher effect size) for each cluster
array = rbind(array,out)
}
}
}
array = array[2:nrow(array),]
array
## Cluster Sample1 Sample2 p-value Cramer's V
## 2 0 Control-P5 Control-P7 0.361085177576349 0.0109394808198687
## 3 1 Control-P5 Control-P7 7.64665654489035e-22 0.191938177732995
## 4 2 Control-P5 Control-P7 5.95614534513094e-29 0.315832817634304
## 5 3 Control-P5 Control-P7 0.000510828865315769 0.100148205635913
## 6 4 Control-P5 Control-P7 0.992014330250716 0.00059079605512425
## 7 5 Control-P5 Control-P7 0.0283761617600838 0.188661309213492
## 8 6 Control-P5 Control-P7 0.457057004447663 0.0968217050180083
## 9 0 Control-P5 Deprived-P5 1.34437556296484e-07 0.0631561954152618
## 10 1 Control-P5 Deprived-P5 0.240570556886315 0.0234524831429395
## 11 2 Control-P5 Deprived-P5 0.0061179336870449 0.0775382950903311
## 12 3 Control-P5 Deprived-P5 0.313212586343204 0.029064601013856
## 13 4 Control-P5 Deprived-P5 5.74325465511809e-61 0.972365568980841
## 14 5 Control-P5 Deprived-P5 0.539737740038465 0.052776681627256
## 15 6 Control-P5 Deprived-P5 0.384777394857507 0.113150763165275
## 16 0 Control-P5 Deprived-P7 0.909627625387819 0.00135959959488482
## 17 1 Control-P5 Deprived-P7 2.18972382631473e-26 0.212405620444991
## 18 2 Control-P5 Deprived-P7 5.61797419449398e-55 0.44170221339753
## 19 3 Control-P5 Deprived-P7 0.923147899499962 0.00278019739973707
## 20 4 Control-P5 Deprived-P7 0.966287651239758 0.00249480889152806
## 21 5 Control-P5 Deprived-P7 3.36736328288764e-06 0.399953288782341
## 22 6 Control-P5 Deprived-P7 0.254690956627095 0.148288976969146
## 23 0 Control-P7 Deprived-P5 5.46787971596413e-11 0.0785468668922356
## 24 1 Control-P7 Deprived-P5 2.58495939617642e-15 0.158062913908414
## 25 2 Control-P7 Deprived-P5 1.04460165868827e-14 0.21874223664928
## 26 3 Control-P7 Deprived-P5 0.0289975590096986 0.0629279919150176
## 27 4 Control-P7 Deprived-P5 7.72486774798787e-83 1.13813382192853
## 28 5 Control-P7 Deprived-P5 0.192595588304233 0.112139134298106
## 29 6 Control-P7 Deprived-P5 0.0769274145317663 0.230278022713877
## 30 0 Control-P7 Deprived-P7 0.297711436963763 0.0124733142024764
## 31 1 Control-P7 Deprived-P7 0.0605471021946722 0.0375057751916102
## 32 2 Control-P7 Deprived-P7 1.22108212629581e-11 0.191702166205112
## 33 3 Control-P7 Deprived-P7 0.001169424801631 0.0935549222150577
## 34 4 Control-P7 Deprived-P7 0.999999999999931 5.08903694414157e-15
## 35 5 Control-P7 Deprived-P7 0.00390611973253633 0.248356813635614
## 36 6 Control-P7 Deprived-P7 0.71897617142499 0.0468455443258497
## 37 0 Deprived-P5 Deprived-P7 4.77246258376725e-07 0.0603120282246738
## 38 1 Deprived-P5 Deprived-P7 1.9647586227062e-19 0.180160327635645
## 39 2 Deprived-P5 Deprived-P7 5.87760093287866e-38 0.364284993161216
## 40 3 Deprived-P5 Deprived-P7 0.393022579752167 0.0246161575174639
## 41 4 Deprived-P5 Deprived-P7 2.43986306811705e-57 0.942070014481652
## 42 5 Deprived-P5 Deprived-P7 0.000197431955903215 0.320362849330336
## 43 6 Deprived-P5 Deprived-P7 0.0368052047501156 0.271823810217491
write.csv(array,file = file.path(dir,"Spreadsheets",paste0(name,"_contingency.csv")),row.names = T,quote = F)